Antes de comenzar, se comprueba que las librerías a emplear a lo largo de los ejercicios estén instaladas y, en caso contrario, proceder a instalarlas con el siguiente código:
comprobar <- function(pkg){
new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if (length(new.pkg))
install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}
paquetes<-c("tidyverse","factoextra","FactoMineR","cluster","readxl",
"ggplot2","corrplot","psych","heatmaply","VIM",
"NbClust","scales","descr","caret","magrittr", "dendextend",
"ggpubr")
comprobar(paquetes)
## tidyverse factoextra FactoMineR cluster readxl ggplot2 corrplot
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## psych heatmaply VIM NbClust scales descr caret
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## magrittr dendextend ggpubr
## TRUE TRUE TRUE
A continuación, antes de proceder a la resolución de las preguntas, se cargarán los datos de las provincias en un dataframe:
# En un script interactivo se podría hacer con el comando choose.files:
# data<-read_xlsx(choose.files(), sheet = "Hoja1")
data<-read_xlsx("Provincias.xlsx", sheet = "Hoja1")
data<-as.data.frame(data)
head(data)
## Prov Poblacion Mortalidad Natalidad IPC NumEmpresas Industria
## 1 Albacete 396987 9.41 9.02 101.819 26701 3332
## 2 Alicante 1868438 8.00 8.52 102.309 130438 9992
## 3 Almería 701688 6.91 11.41 101.454 40327 2191
## 4 Álava 321932 7.72 10.26 102.743 19548 2048
## 5 Asturias 1061756 12.16 6.26 102.070 67451 3496
## 6 Badajoz 690929 9.54 8.96 101.036 39640 3059
## Construccion CTH Infor AFS APT TasaActividad TasaParo Ocupados PIB
## 1 3428 11111 189 624 3278 57.72 23.28 144.8 7235991
## 2 17418 52720 1851 2745 19865 58.21 23.39 687.1 32139347
## 3 5388 18184 366 944 5840 60.63 31.08 234.4 11909684
## 4 2669 7557 323 354 3341 57.86 13.33 133.0 10716021
## 5 8435 28045 787 1445 10598 51.42 16.97 389.0 21770433
## 6 4599 18420 332 948 5235 56.00 29.63 224.6 10581366
## CANE TVF VS
## 1 20888 215948 30244
## 2 25968 1274096 326705
## 3 22945 395086 72486
## 4 3691 155767 9791
## 5 23910 613905 73250
## 6 37438 372493 49441
Se ha comprobado que los datos cargados son un dataframe por lo que se evaluará su estructura interna para hacerse una idea de su disposición y de los datos faltantes:
summary(data)
## Prov Poblacion Mortalidad Natalidad
## Length:52 Min. : 84509 Min. : 5.820 Min. : 5.550
## Class :character 1st Qu.: 322203 1st Qu.: 7.855 1st Qu.: 7.700
## Mode :character Median : 614723 Median : 9.120 Median : 8.975
## Mean : 899449 Mean : 9.379 Mean : 8.839
## 3rd Qu.:1019030 3rd Qu.:10.688 3rd Qu.: 9.610
## Max. :6454440 Max. :14.360 Max. :19.330
## IPC NumEmpresas Industria Construccion
## Min. :100.6 Min. : 3749 Min. : 75 Min. : 309
## 1st Qu.:101.9 1st Qu.: 22822 1st Qu.: 1704 1st Qu.: 2972
## Median :102.3 Median : 38000 Median : 2516 Median : 5070
## Mean :102.4 Mean : 61286 Mean : 3808 Mean : 7805
## 3rd Qu.:102.7 3rd Qu.: 65083 3rd Qu.: 4106 3rd Qu.: 8350
## Max. :104.8 Max. :508612 Max. :27416 Max. :59661
## CTH Infor AFS APT
## Min. : 2030 Min. : 35.0 Min. : 50.0 Min. : 504
## 1st Qu.: 9243 1st Qu.: 185.5 1st Qu.: 486.8 1st Qu.: 3091
## Median : 15488 Median : 369.5 Median : 813.5 Median : 5440
## Mean : 23741 Mean : 1131.9 Mean : 1378.3 Mean : 10854
## 3rd Qu.: 27567 3rd Qu.: 868.2 3rd Qu.: 1429.2 3rd Qu.: 10627
## Max. :158331 Max. :19058.0 Max. :12357.0 Max. :123863
## TasaActividad TasaParo Ocupados PIB
## Min. :47.41 Min. :11.95 Min. : 24.6 Min. : 1397441
## 1st Qu.:55.45 1st Qu.:15.39 1st Qu.: 132.4 1st Qu.: 6509393
## Median :57.79 Median :19.51 Median : 222.4 Median : 11883640
## Mean :57.84 Mean :21.17 Mean : 347.1 Mean : 20275145
## 3rd Qu.:60.07 3rd Qu.:27.75 3rd Qu.: 389.8 3rd Qu.: 21242697
## Max. :68.69 Max. :37.18 Max. :2806.4 Max. :198652445
## CANE TVF VS
## Min. : 3 Min. : 26233 Min. : 200
## 1st Qu.: 9758 1st Qu.: 211628 1st Qu.: 38941
## Median :14037 Median : 335934 Median : 56412
## Mean :19035 Mean : 484781 Mean : 70799
## 3rd Qu.:26020 3rd Qu.: 532066 3rd Qu.: 80625
## Max. :68037 Max. :2894679 Max. :326705
aggr_plot <- aggr(data, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(data), cex.axis=.5, gap=1, ylab=c("Histogram of missing data","Pattern"))
##
## Variables sorted by number of missings:
## Variable Count
## Prov 0
## Poblacion 0
## Mortalidad 0
## Natalidad 0
## IPC 0
## NumEmpresas 0
## Industria 0
## Construccion 0
## CTH 0
## Infor 0
## AFS 0
## APT 0
## TasaActividad 0
## TasaParo 0
## Ocupados 0
## PIB 0
## CANE 0
## TVF 0
## VS 0
El gráfico muestra que en el dataset no hay datos faltantes, por lo que no será necesario proceder a su inserción.
Por comodidad, se añadirán los nombres de las provincias de la col.1 como rownames eliminando esta columna del conjunto tras la operación.
rownames(data)<-data$Prov
data<-data[,-1]
Para responder a esta pregunta, se creará la matriz de correlaciones y se representará mediante el paquete corrplot:
matriz_cor<-cor(data, method = "pearson")
corrplot(matriz_cor, type="upper", order="hclust",
tl.col="black", tl.cex=0.6, tl.srt=90)
Con este gráfico se puede apreciar que las variables más correlacionadas negativamente son Mortalidad con TasaActividad, Natalidad y TasaParo así como la TasaParo con el IPC. Una interpretación de esto es que en las provincias con mayor mortalidad, la natalidad suele ser baja como pudiera ser en las regiones más envejecidas donde hay fallecimientos, pero pocas familias jóvenes.
Para esta pregunta se utilizará la función PCA de factominer combinada después con fviz de factoextra:
pca<-PCA(data,ncp = 7, scale.unit = TRUE, graph = TRUE)
# Especificando el argumento scale.unit=TRUE el PCA se realiza sobre las variables estandarizadas que es como aplicarlo sobre la matriz de correlaciones
eig<-get_eigenvalue(pca)
Observando los 2 primeros gráficos se extrae que la componente 1 explica el 63,7% de la variabilidad de los datos y la componente 2 el 14,23%, algo muy positivo porque gran parte de la muestra podría quedar reducida a estas 2.
Por otro lado, en cuanto a las provincias como tal, Madrid y Barcelona están muy alejadas del resto de provincias y muy ligadas a la PC1, mientras que hay luego otros 2 grupos más explicados por la PC2, uno con provincias del sur como Melilla, Málaga, Sevilla… y otro con más del centro-norte como Palencia, Pontevedra, Álava, Zamora, etc.
Comprobando los valores de los autovalores:
knitr:: kable(eig, digits = 2, caption="Autovalores")
| eigenvalue | variance.percent | cumulative.variance.percent | |
|---|---|---|---|
| Dim.1 | 11.47 | 63.70 | 63.70 |
| Dim.2 | 2.56 | 14.23 | 77.93 |
| Dim.3 | 1.63 | 9.08 | 87.01 |
| Dim.4 | 0.93 | 5.19 | 92.19 |
| Dim.5 | 0.46 | 2.54 | 94.73 |
| Dim.6 | 0.41 | 2.30 | 97.03 |
| Dim.7 | 0.31 | 1.71 | 98.74 |
| Dim.8 | 0.12 | 0.65 | 99.39 |
| Dim.9 | 0.07 | 0.41 | 99.79 |
| Dim.10 | 0.02 | 0.11 | 99.91 |
| Dim.11 | 0.01 | 0.05 | 99.96 |
| Dim.12 | 0.00 | 0.02 | 99.98 |
| Dim.13 | 0.00 | 0.01 | 99.99 |
| Dim.14 | 0.00 | 0.00 | 99.99 |
| Dim.15 | 0.00 | 0.00 | 100.00 |
| Dim.16 | 0.00 | 0.00 | 100.00 |
| Dim.17 | 0.00 | 0.00 | 100.00 |
| Dim.18 | 0.00 | 0.00 | 100.00 |
Atendiendo a la tabla, con las 4 primeras componentes, se podría explicar el 92,19% de los datos siendo que el valor de la 4 está muy cerca de 1. El número más adecuado de componentes se puede extraer con el scree plot:
fviz_eig(pca, addlabels=TRUE)
El número más óptimo serían 3 o 4 componentes. Como se ha explicado antes, se tomarán 4 para superar el 90% de datos explicados.
En el apartado anterior se concluyó que 4 componentes serían suficientes de manera que:
pca_4<-PCA(data,ncp = 4, scale.unit = TRUE, graph = TRUE)
knitr::kable(pca_4$svd$V, digits=3, caption = "Autovectores")
| 0.294 | 0.002 | 0.050 | -0.053 |
| -0.106 | -0.527 | 0.189 | -0.161 |
| 0.041 | 0.495 | -0.271 | -0.110 |
| 0.110 | -0.365 | -0.262 | 0.435 |
| 0.294 | -0.026 | 0.008 | -0.069 |
| 0.286 | -0.045 | 0.046 | 0.023 |
| 0.293 | -0.045 | -0.012 | -0.026 |
| 0.293 | -0.011 | 0.049 | -0.028 |
| 0.282 | -0.042 | -0.065 | -0.222 |
| 0.292 | -0.016 | 0.040 | -0.092 |
| 0.291 | -0.029 | -0.028 | -0.142 |
| 0.114 | 0.331 | -0.363 | 0.463 |
| -0.014 | 0.462 | 0.387 | -0.220 |
| 0.294 | -0.017 | 0.002 | -0.060 |
| 0.291 | -0.036 | -0.037 | -0.134 |
| 0.018 | 0.096 | 0.657 | 0.278 |
| 0.292 | -0.002 | 0.100 | 0.044 |
| 0.172 | 0.048 | 0.290 | 0.567 |
Cada una de las 4 columnas mostradas representa una de las 4 componenetes escogidas para el análisis y las líneas son las variables originales del dataset con lo que la expresión de la PC1 sería:
$ PC1 = 0,294Poblacion - 0,106Mortalidad + 0,041Natalidad + … + 0,292TVF + 0,172VS $
var<-get_pca_var(pca_4)
knitr::kable(var$cor, digits=2, caption = "Correlaciones PC con variables")
| Dim.1 | Dim.2 | Dim.3 | Dim.4 | |
|---|---|---|---|---|
| Poblacion | 0.99 | 0.00 | 0.06 | -0.05 |
| Mortalidad | -0.36 | -0.84 | 0.24 | -0.16 |
| Natalidad | 0.14 | 0.79 | -0.35 | -0.11 |
| IPC | 0.37 | -0.58 | -0.34 | 0.42 |
| NumEmpresas | 1.00 | -0.04 | 0.01 | -0.07 |
| Industria | 0.97 | -0.07 | 0.06 | 0.02 |
| Construccion | 0.99 | -0.07 | -0.02 | -0.03 |
| CTH | 0.99 | -0.02 | 0.06 | -0.03 |
| Infor | 0.95 | -0.07 | -0.08 | -0.21 |
| AFS | 0.99 | -0.03 | 0.05 | -0.09 |
| APT | 0.98 | -0.05 | -0.04 | -0.14 |
| TasaActividad | 0.39 | 0.53 | -0.46 | 0.45 |
| TasaParo | -0.05 | 0.74 | 0.50 | -0.21 |
| Ocupados | 1.00 | -0.03 | 0.00 | -0.06 |
| PIB | 0.99 | -0.06 | -0.05 | -0.13 |
| CANE | 0.06 | 0.15 | 0.84 | 0.27 |
| TVF | 0.99 | 0.00 | 0.13 | 0.04 |
| VS | 0.58 | 0.08 | 0.37 | 0.55 |
De cara a ver las correlaciones más fuertes de las variables con las componentes, se ordenará el df de var$cor de mayor a menor para cada componente y se extraerán las variables con mayor correlación (positiva o negativa):
df<-as.data.frame(var$cor)
df %>% select("Dim.1") %>% arrange(desc(Dim.1))
## Dim.1
## Ocupados 0.99699549
## NumEmpresas 0.99616251
## Poblacion 0.99394316
## Construccion 0.99305566
## CTH 0.99170521
## AFS 0.99035020
## TVF 0.98730511
## PIB 0.98507163
## APT 0.98445474
## Industria 0.96717807
## Infor 0.95322382
## VS 0.58359002
## TasaActividad 0.38716977
## IPC 0.37218373
## Natalidad 0.13757149
## CANE 0.06023588
## TasaParo -0.04739450
## Mortalidad -0.35992669
En la PC1 tienen una correlación casi perfecta los Ocupados, NumEmpresas, Poblaciónm Construcción, etc. hasta Infor teniendo las correlacionadas negativamente valores muy bajos
df %>% select(Dim.2) %>% arrange(desc(Dim.2))
## Dim.2
## Natalidad 0.792749686
## TasaParo 0.738572741
## TasaActividad 0.529035826
## CANE 0.153910432
## VS 0.076521234
## Poblacion 0.003989744
## TVF -0.002559620
## CTH -0.016989690
## AFS -0.026329340
## Ocupados -0.026944763
## NumEmpresas -0.041840166
## APT -0.047097804
## PIB -0.057858046
## Infor -0.067235670
## Industria -0.071587054
## Construccion -0.072618036
## IPC -0.584573656
## Mortalidad -0.843475898
A la PC2 contribuyen positivamente Natalidad, TasaParo y TasaActividad mientras que de forma negativa, Mortalidad e IPC son las que más destacan.
df %>% select(Dim.3) %>% arrange(desc(Dim.3))
## Dim.3
## CANE 0.839526520
## TasaParo 0.495165389
## VS 0.370894920
## Mortalidad 0.241432113
## TVF 0.127527357
## Poblacion 0.063848580
## CTH 0.062434517
## Industria 0.059274514
## AFS 0.050562974
## NumEmpresas 0.010165492
## Ocupados 0.003105467
## Construccion -0.015444759
## APT -0.035880121
## PIB -0.047926363
## Infor -0.082630799
## IPC -0.335325313
## Natalidad -0.346016593
## TasaActividad -0.463689989
A la PC3 CANE es la mejor contribuyente en la parte positiva y TasaActividad negativamente aunque sin llegar al 0,5.
df %>% select(Dim.4) %>% arrange(desc(Dim.4))
## Dim.4
## VS 0.54810069
## TasaActividad 0.44752178
## IPC 0.42045677
## CANE 0.26905683
## TVF 0.04266413
## Industria 0.02266458
## Construccion -0.02550892
## CTH -0.02665589
## Poblacion -0.05160463
## Ocupados -0.05803848
## NumEmpresas -0.06681923
## AFS -0.08860621
## Natalidad -0.10599223
## PIB -0.12912629
## APT -0.13710328
## Mortalidad -0.15590654
## TasaParo -0.21281673
## Infor -0.21409732
Para la PC4, VS es la única destacable.
En la primera visualización se comprobará cuál es la correlación respecto a las 2 primeras componentes:
fviz_pca_var(pca_4,
axes = c(1,2),
col.var = "cos2", # Color by contributions to the PC
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping
)
La mayoría de las variables están muy fuertemente correlacionadas de forma positiva con la componente 1 como son Población, el PIB, Industria, Construcción, etc. Mientras, en la componente 2, la Mortalidad posee una alta correlación negativa con la misma, a diferencia de la Natalidad o la TasaParo que se correlacionan de forma positiva. Por tanto, la componente 1 recoge variables que muestran el número de empresas de diferentes sectores (variables más puramente económicas) a la par que la componente 2 explica variables sociales como Natalidad, Mortalidad, tasas de paro, etc.
En cuanto a las componentes 3 y 4:
fviz_pca_var(pca_4,
axes = c(3,4),
col.var = "cos2", # Color by contributions to the PC
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping
)
En este gráfico hay bastante menos correlación que en el anterior destacando en la componente 3 la representación de indicadores socioeconómicos como el paro, la actividad o el IPC, así como las explotaciones agrarias. En la componente 4 no hay apenas representación aunque podría solaparse con la interpretación de la 3 al formar las flechas casi 45º entre componentes.
variables es la que está peor explicada?
Esta pregunta se puede responder comporbando los cos2 de cada variable pues muestra la proporción de la varianza de cada variable que es explicada por cada componente:
knitr::kable(var$cos2, digits = 2, caption="Cosenos al cuadrado")
| Dim.1 | Dim.2 | Dim.3 | Dim.4 | |
|---|---|---|---|---|
| Poblacion | 0.99 | 0.00 | 0.00 | 0.00 |
| Mortalidad | 0.13 | 0.71 | 0.06 | 0.02 |
| Natalidad | 0.02 | 0.63 | 0.12 | 0.01 |
| IPC | 0.14 | 0.34 | 0.11 | 0.18 |
| NumEmpresas | 0.99 | 0.00 | 0.00 | 0.00 |
| Industria | 0.94 | 0.01 | 0.00 | 0.00 |
| Construccion | 0.99 | 0.01 | 0.00 | 0.00 |
| CTH | 0.98 | 0.00 | 0.00 | 0.00 |
| Infor | 0.91 | 0.00 | 0.01 | 0.05 |
| AFS | 0.98 | 0.00 | 0.00 | 0.01 |
| APT | 0.97 | 0.00 | 0.00 | 0.02 |
| TasaActividad | 0.15 | 0.28 | 0.22 | 0.20 |
| TasaParo | 0.00 | 0.55 | 0.25 | 0.05 |
| Ocupados | 0.99 | 0.00 | 0.00 | 0.00 |
| PIB | 0.97 | 0.00 | 0.00 | 0.02 |
| CANE | 0.00 | 0.02 | 0.70 | 0.07 |
| TVF | 0.97 | 0.00 | 0.02 | 0.00 |
| VS | 0.34 | 0.01 | 0.14 | 0.30 |
corrplot(var$cos2, is.corr = FALSE, tl.cex=0.6, tl.col="black",
cl.ratio = 1)
La componente 1 explica la práctica totalidad de múltiples variables mencionadas en el apartado previo (Industria, Infor, PIB, Ocupados, etc.) y la 2 aquellas no tan explicadas por la 1, de manera que se compensan esos nichos y por ello el poder explicativo es tan alto.
fviz_cos2(pca_4, choice = "var", axes=1:4, tl.cex=0.6)
La variable peor explicada es el IPC, con una varianza total explicada por las 4 componentes de \(0.14+0.34+0.11+0.18 = 0.77\)
variables es la que está peor explicada?
En este caso se deberá atender al apartado contrib (contribuciones) del objeto var:
knitr::kable(var$contrib, digits = 2, caption="Contribuciones de las variables")
| Dim.1 | Dim.2 | Dim.3 | Dim.4 | |
|---|---|---|---|---|
| Poblacion | 8.62 | 0.00 | 0.25 | 0.29 |
| Mortalidad | 1.13 | 27.79 | 3.57 | 2.60 |
| Natalidad | 0.17 | 24.54 | 7.33 | 1.20 |
| IPC | 1.21 | 13.35 | 6.88 | 18.93 |
| NumEmpresas | 8.65 | 0.07 | 0.01 | 0.48 |
| Industria | 8.16 | 0.20 | 0.22 | 0.05 |
| Construccion | 8.60 | 0.21 | 0.01 | 0.07 |
| CTH | 8.58 | 0.01 | 0.24 | 0.08 |
| Infor | 7.92 | 0.18 | 0.42 | 4.91 |
| AFS | 8.55 | 0.03 | 0.16 | 0.84 |
| APT | 8.45 | 0.09 | 0.08 | 2.01 |
| TasaActividad | 1.31 | 10.93 | 13.16 | 21.44 |
| TasaParo | 0.02 | 21.30 | 15.00 | 4.85 |
| Ocupados | 8.67 | 0.03 | 0.00 | 0.36 |
| PIB | 8.46 | 0.13 | 0.14 | 1.79 |
| CANE | 0.03 | 0.93 | 43.13 | 7.75 |
| TVF | 8.50 | 0.00 | 1.00 | 0.19 |
| VS | 2.97 | 0.23 | 8.42 | 32.16 |
corrplot(var$contrib, is.corr = FALSE, tl.cex=0.6, tl.col="black",
cl.ratio = 1)
Las contribuciones de las variables a la variabilidad de las componentes es algo más escueta que viceversa, destacando fundamentalmente las aportaciones de CANE a la componente 3, mortalidad y natalidad a la 2 y VS a la 4 de forma individual.
fviz_contrib(pca_4, choice = "var", axes=1, tl.cex=0.6)
En la componente 1, ninguna variable explica más del 10% de su variabilidad, siendo la que menos aporta la Tasaparo, confirmando que en esta componente contribuyen más las variables puramente económicas que las socioeconómicas o demográficas.
fviz_contrib(pca_4, choice = "var", axes=2, tl.cex=0.6)
A la componente 2, por su carácter más demógrafico y social, ortogonal a la 1, las que menos aportan son variables como TVF o la Población
fviz_contrib(pca_4, choice = "var", axes=3, tl.cex=0.6)
A la componente 3 contribuyen de forma nula los Ocupados o NumEmpresas
fviz_contrib(pca_4, choice = "var", axes=4, tl.cex=0.6)
En la componente 4 las variables que menos aportan con Construcción e Industria
socioeconómicos para estas provincias?
Para responder a ello se presentarán una serie de biplots con los individuos (las provincias) sobre los ejes de las componentes de la 1 a la 4:
fviz_pca_ind(pca_4, col.ind = "cos2",
axes = c(1,2),
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
col.cex=0.2,
repel = TRUE # Avoid text overlapping (slow if many points)
) +xlim(-10,20)
En este gráfico se observa que en la PC1, Barcelona y Madrid tienen un valor muy superior a todas las demás provincias por ser los motores económicos y destacar en altos índices de empresas de todo tipo y población. Todas las demás se encuentran muy cerca del origen de coordenadas por lo que ninguna es excesivamente baja en estos aspectos económicos salvo aquellas más rurales, que conforman la España vaciada como Soria, Teruel, Albacete, etc.
En la PC2, la componente de métricas más sociales, destacan com valores positivos y altos Ceuta y Melilla, Cádiz, Almería, Málaga… en general provincias más del sur que presentan unas tasas de natalidad y paro y superiores a zonas del centro y norte como Zamora, Ourense, Lugo, Cantabria, etc. que poseen una mayor mortalidad por una población más envejecida, pero con menor nivel de desempleo.
Para las componentes 3 y 4:
fviz_pca_ind(pca_4, col.ind = "cos2",
axes = c(3,4),
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
col.cex=0.2,
repel = TRUE # Avoid text overlapping (slow if many points)
) +xlim(-10,20)
Madrid y Barcelona ya no se encuentran completamente separadas en la componente horizontal (en este caso la 3), ni tampoco ninguna otra provincia se aleja del origen de coordenadas, mostrando poca relación con las explotaciones agrarias o una tasa de paro fuera de lo normal. En cuanto a la componente 4, ciudades como Madrid, Melilla y Ceuta se encuentran a parte por una mayor natalidad o tasa de paro.
A continuación, se representarán conjuntamente las observaciones con las variables dispuestas en los ejes:
fviz_pca_biplot(pca_4, repel=TRUE, col.var="#2E9FD9",
col.cex=0.2, col.ind = "#696969")+xlim(-10,20)
En esta representación se puede ver provincias como Lugo, Zamora, Ourense o Teruel tienen una alta mortalidad. Melilla, Ceuta y Cádiz destacan por su natalidad. Sevilla, Murcia y Ávila tienen buena tasa de actividad y las ya mencionadas Madrid, Barcelona y en menor medida Valencia, Alicante y Navarra son los motores económicos.
fviz_pca_biplot(pca_4, repel=TRUE, axes=c(3,4), col.var="#2E9FD9",
col.cex=0.2, col.ind = "#696969")+xlim(-10,20)
desarrollo económico de una provincia, como se podría construir utilizando una combinación lineal de todas las variables. ¿Cuál sería el valor de dicho índice en Madrid? ¿Cual sería su valor en Melilla?
En los anteriores apartados se ha comprobado que la componente 1 es la que reúne las variables económicas, por lo que esta podría ser un buen índice resultado de la combinación lineal de las variables de forma estandarizada. Es decir, la PC1 viene dada por la expresión: $ PC1 = 0,294Poblacion - 0,106Mortalidad + 0,041Natalidad + … + 0,292TVF + 0,172VS $ en la que cada para cada provincia se tomará el valor que esta presenta en alguna de las variables y se le restará la media para el valor de dicha variable, todo ello dividido entre la desviación típica de la variable. En el caso de Madrid por ejemplo hay una población de 6454440 a la que se restará la media de 899448 y se dividirá por 1154297:
mean(data$Poblacion, na.rm = TRUE)
## [1] 899448.9
sd(data$Poblacion, na.rm = TRUE)
## [1] 1154297
Después se multiplicará el resultado por 0,294. Esto se repetiría para cada variable y cada provincia, pero se puede extraer directamente de:
ind<-get_pca_ind(pca_4)
knitr::kable(ind$coord, digits =3,caption = "Valores de los individuo
s en las Cp")
| Dim.1 | Dim.2 | Dim.3 | Dim.4 | |
|---|---|---|---|---|
| Albacete | -1.410 | 0.473 | 0.096 | -0.458 |
| Alicante | 3.384 | 0.540 | 1.919 | 2.420 |
| Almería | -0.617 | 2.614 | 0.208 | -0.142 |
| Álava | -1.444 | -0.001 | -2.032 | -0.113 |
| Asturias | -0.204 | -1.953 | 1.298 | -0.714 |
| Badajoz | -1.048 | 1.168 | 1.796 | -0.854 |
| Balears | 1.526 | 0.260 | -2.519 | 2.364 |
| Barcelona | 13.683 | -1.612 | -0.867 | -0.279 |
| Bizkaia | 0.576 | -1.508 | -1.180 | -0.383 |
| Burgos | -1.202 | -1.550 | -0.892 | 0.989 |
| Cantabria | -0.849 | -1.126 | -0.594 | 0.401 |
| Castellón | -0.690 | 0.821 | 0.518 | 0.382 |
| Ceuta | -2.125 | 3.326 | -1.811 | -1.165 |
| Ciudad Real | -1.392 | 0.815 | 1.819 | -0.774 |
| Coruña | 0.635 | -1.442 | 0.759 | 0.244 |
| Cuenca | -2.132 | -0.569 | 1.048 | -0.932 |
| Cáceres | -1.503 | -0.180 | 1.269 | -0.223 |
| Cádiz | 0.128 | 1.771 | 0.369 | -0.249 |
| Córdoba | -0.594 | 0.811 | 1.292 | -0.169 |
| Gipuzkoa | -0.281 | -1.485 | -1.879 | 0.266 |
| Girona | 0.388 | 0.544 | -1.387 | 1.638 |
| Granada | -0.072 | 1.124 | 1.511 | 0.500 |
| Guadalajara | -1.408 | 1.542 | -1.849 | 0.906 |
| Huelva | -1.265 | 1.168 | 0.013 | -0.327 |
| Huesca | -1.776 | -0.984 | -0.587 | 0.304 |
| Jaén | -1.221 | 1.424 | 3.407 | -0.479 |
| León | -1.464 | -2.016 | 0.877 | -0.809 |
| Lleida | -0.835 | -0.136 | -1.472 | 1.185 |
| Lugo | -1.826 | -2.785 | 0.871 | -0.353 |
| Madrid | 16.778 | -0.366 | -0.849 | -2.365 |
| Melilla | -2.218 | 4.782 | -1.905 | -2.231 |
| Murcia | 1.522 | 1.442 | 0.580 | 0.901 |
| Málaga | 2.006 | 1.325 | 0.869 | 1.104 |
| Navarra | -0.653 | 0.078 | -1.096 | -0.081 |
| Ourense | -1.965 | -2.858 | 1.098 | -1.204 |
| Palencia | -2.122 | -1.951 | -0.695 | -0.142 |
| Palmas | 0.092 | 1.857 | -0.330 | -0.824 |
| Pontevedra | 0.036 | -0.607 | 0.052 | -0.328 |
| Rioja | -1.383 | -0.484 | -1.354 | 0.288 |
| Salamanca | -1.612 | -1.425 | -0.121 | -0.086 |
| Santa Cruz | -0.029 | 1.573 | -0.126 | -0.432 |
| Segovia | -1.931 | -1.015 | -1.110 | 0.263 |
| Sevilla | 1.948 | 1.775 | 0.712 | -0.546 |
| Soria | -2.399 | -1.857 | -0.778 | -0.503 |
| Tarragona | 0.175 | 1.040 | -0.102 | 1.512 |
| Teruel | -2.185 | -1.082 | -0.351 | -0.413 |
| Toledo | -0.461 | 1.449 | 0.812 | 0.463 |
| Valencia | 4.770 | 0.360 | 2.961 | 2.349 |
| Valladolid | -1.007 | -0.704 | -1.052 | 0.196 |
| Zamora | -2.287 | -3.169 | 0.578 | -0.622 |
| Zaragoza | 0.115 | -0.237 | -0.156 | 0.068 |
| Ávila | -2.152 | -0.982 | 0.365 | -0.545 |
Madrid tendría un valor respecto a potencia económica de 16,778 mientras que Melilla presentaría un -2,218 haciendo notar numéricamente la diferencia abismal entre ambas que quedó manifiesta en los gráficos.
Se comienza creando una nueva versión de los datos (data_st) estandarizados a partir del dataset original y eliminando las variables no numéricas:
data_st<-scale(data)
A continuación se crean mapas de calor para ambos conjuntos:
heatmaply(data, seriate="mean",
row_dend_left = TRUE, plot_method = "plotly")
heatmaply(data_st, seriate="mean",
row_dend_left = TRUE, plot_method = "plotly")
Con los datos no estandarizados, parece muy clara la división en 3 grupos del conjunto de datos, Madrid y barcelona por un lado, Las potencias turísticas e industriales del mediterráneo y cantábrico por otro y todas las demás en un tercer grupo (aunque podría hacerse una cuarta segregación dando lugar a 4 grupos). No obstante, al no estandarizar y permitir que las variables con distintas unidades de medida se mezclen (población total y tasa de actividad en % por ejemplo) los resultados no reflejan bien las relaciones.
En el segundo gráfico con estas ya estandarizadas aparece una distinción entre grupos menos clara, lo que puede lleva a deducir que hay más grupos: Madrid y Barcelona, las ciudades autónomas, Valencia y Alicante y luegos 2 grupos adicionales con el resto de comunidades.
Para llevar a cabo el análisis cluster, una vez que se han explorado superficialmente los datos con el mapa de calor, se debe crear la matriz de distancias entre las variables (únicamente sobre ellas estandarizadas pues es la manera correcta de hacerlo):
dist_st <- dist(data_st, method = "euclidean")
Con el anterior comando se extrae la distancia euclídea.
A continuación se procede a su visualización:
fviz_dist(dist_st)
También se puede crear la visualización mediante el heatmap visto previamente y agrupando los más parecidos.
heatmaply(as.matrix(dist_st), seriate="mean",
row_dend_left = TRUE, plot_method = "plotly")
Las provincias más agrupadas según se refleja en estos gráficos son Madrid y Barcelona, muy cerca del origen de coordenadas, luego se distingue un gran grupo que incuye provincias como Huesca, Teruel, Zamora, Lugo, Bizkaia, Alava, etc. Un pequeño grupo con Guadalajara, Girona, Baleares y Tarragona. Ceuta y Melilla en su grupo particular y luego algunas otras.
Por último, se procede a representar el dendrograma de los datos estandarizados mediante el método ward.D2:
res_hc <- hclust(dist_st, method="ward.D2")
fviz_dend(res_hc, cex = 0.5)
Según el dendrograma, podría haber distinto número de clusters. La opción más adecuada sin sobreestimar es 3 clusters, a la altura de entre 10 y 15, o entrando en cada uno de los 2 grupos más a la derecha podrían darse hasta 6 como se muestra en las siguientes representaciones:
# 3 clusters
as.dendrogram(res_hc) %>% set("branches_k_color",
value = c("red", "blue", "green"), k = 3) %>%
plot(main= "3 clusters")
# 6 clusters
as.dendrogram(res_hc) %>% set("branches_k_color",
value = c("red", "blue", "green", "brown", "orange","purple"), k = 6) %>%
plot(main= "3 clusters")
Se deberá comprobar con los métodos de optimización cuál es el mejor número.
La opción más razonable es hacer 3 clusters por lo que se tomará este como número óptimo por el momento. Se comprueba cuántos individuos hay por cluster:
group<-cutree(res_hc, k=3)
knitr::kable(table(group), caption="Nº individuos por cluster")
| group | Freq |
|---|---|
| 1 | 31 |
| 2 | 19 |
| 3 | 2 |
Se visualizan estos grupos representados en el dendrograma:
fviz_dend(res_hc, k=3, cex=0.5, color_labels_by_k = TRUE, rect=TRUE)
También se pueden visualizar los cluster sobre las componentes principales:
fviz_cluster(list(data = data_st,cluster=group),
ellipse.type = "convex",
repel = TRUE,
ggtheme = theme_minimal(),
show.clust.cent=FALSE,
main = "Representación sobre PC1 y PC2"
)
El número óptimo de clusters según el criterio Elbow y de Silhouette sería:
# Elbow method
fviz_nbclust(data_st, kmeans, method = "wss") +
geom_vline(xintercept = c(3,6), linetype=2)
labs(subtitle = "Elbow method")
## $subtitle
## [1] "Elbow method"
##
## attr(,"class")
## [1] "labels"
# Silhouette method
fviz_nbclust(data_st, kmeans, method = "silhouette")+
labs(subtitle = "Silhouette method")
El criterio de Elbow muestra 2 puntos de inflexión claros en n=3 y n=6, tal como se mencionaba en apartados anteriores. En principio 3 sigue siendo más estable, pero se deberá probar a visualizar con n=6 por si arroja resultados más adecuados.
Con el criterio de Silhouette ocurre algo similar sólo que los picos se encuentran fundamentalmente en 3 y 5.
El agrupamiento no jerárquico se llevará a cabo con el algoritmo de kmeans y, antes de proseguir, se comprobará con silhouette la calidad con 3 y 6 clusters:
set.seed(123)
km_res_3<-kmeans(data_st, 3)
km_res_6<-kmeans(data_st, 6)
sil_3 <- silhouette(km_res_3$cluster, dist(data_st))
rownames(sil_3) <- rownames(data_st)
sil_6 <- silhouette(km_res_6$cluster, dist(data_st))
rownames(sil_6) <- rownames(data_st)
fviz_silhouette(sil_3)
## cluster size ave.sil.width
## 1 1 2 0.66
## 2 2 23 0.08
## 3 3 27 0.41
fviz_silhouette(sil_6)
## cluster size ave.sil.width
## 1 1 2 0.46
## 2 2 16 0.16
## 3 3 13 0.36
## 4 4 2 0.60
## 5 5 4 0.21
## 6 6 15 0.24
Tanto con 3 como con 6 cluster se obtiene un valor de with muy similar, con algunos valores negativos (incorrectamente asignados en ambos). Al ejecutar varias veces (y por tanto variar la selección aleatoria con kmeans) los resultados muestran grupos poco constantes, a veces sin widths negativos, grupos muy reducidos o muy extensos.
Por ello, se decide crear una función que muestre cuál es la media de las widths según 3,5 o 6 clusters tras ejecutar 100 veces el algoritmo:
get_average_width<-function(clusters){
widths=0
counter=0
for (i in c(1:100)) {
res<-kmeans(data_st, clusters)
sil <- silhouette(res$cluster, dist(data_st))
rownames(sil) <- rownames(data_st)
widths=widths+mean(sil[,3])
counter=counter+1
}
average=widths/counter
return(average)
}
get_average_width(3)
## [1] 0.2819402
get_average_width(5)
## [1] 0.2297499
get_average_width(6)
## [1] 0.2301892
Tras esta comprobación adicional, 3 clusters ofrecen la mayor width media, de manera que se tomará este número como óptimo para continuar con el análisis.
Se procede a representarlos sobre las componentes principales:
fviz_cluster(list(data = data_st,cluster=group),
ellipse.type = "convex",
repel = TRUE,
ggtheme = theme_minimal(),
show.clust.cent=FALSE,
main = "Representación sobre PC1 y PC2"
) +xlim(-10,20)
La interpretación de esta representación es la misma que se hacía con el análisis de las componentes. Cuanto más a la derecha en el gráfico, mayor es la concentración empresarial, poblacional y en general, mejor situación económica presentan las provincias. Por otro lado,cuanto más hacia arriba se encuentren las provincias, gozarán de mejor situación social en cuanto a natalidad o tasa de actividad mientras que hacia abajo supone mayor mortalidad.
Sobre PC3 y PC3:
fviz_cluster(list(data = data_st,cluster=group),
axes = c(3,4),
ellipse.type = "convex",
repel = TRUE,
ggtheme = theme_minimal(),
show.clust.cent=FALSE,
main = "Representación sobre PC3 y PC4"
) +xlim(-10,20)
Las componentes PC3 y PC4 eran algo menos específicas que las 2 primeras destacando magnitudes como las explotaciones agrarias, la tasa de actividad y el IPC. Esta agrupación de conjuntos muestra que cuanto más a la derecha en el gráfico, más agricultura hay en las provincias como es en los casos de Alicante, Valencia o Murcia y cuanto más a la izquierda y en la parte positiva del eje Y se ubican empresas con buena tasa de actividad y un IPC moderadamente elevado como Girona, Tarragona o Huelva.
Parte de este apartado se realizaó al comienzo del apartado, pero de manera algo más general para averiguar el número de clusters correcto. A continuación se llevará a cabo de forma más minuciosa sabiendo que se dispone de 3 clusters:
set.seed(123)
sil_3 <- silhouette(km_res_3$cluster, dist(data_st))
rownames(sil_3) <- rownames(data_st)
head(sil_3[,1:3])
## cluster neighbor sil_width
## Albacete 3 2 0.23404690
## Alicante 2 3 0.12861643
## Almería 2 3 0.23985301
## Álava 3 2 0.33700472
## Asturias 3 2 0.38504449
## Badajoz 2 3 0.03338686
fviz_silhouette(sil_3)
## cluster size ave.sil.width
## 1 1 2 0.66
## 2 2 23 0.08
## 3 3 27 0.41
Con el primer head se muestra el valor del width que posee cada una de las provincias siendo que cuanto más cerca de 1 esté, mejor agrupado estará el individuo y si está próximo a 0, podría ubicarse en cualquier otro grupo. Con la función para representar el gráfico se muestra la media de esta medida para cada uno de los clusters con sus correspondientes observaciones. Los cluster de 27 y 2 individuos tienen un valor muy elevado, pero el de 23 miembros se ve mermado por las provincias con valor negativo (Girona, Huelva, Guadalajara, Baleares y Castellón y Ciudad Real) que, pese a todo, están muy cerca de 0.
sil_3[,1:3] %>% as.data.frame() %>% filter(sil_width<0)
## cluster neighbor sil_width
## Balears 2 3 -0.03186983
## Castellón 2 3 -0.06554316
## Ciudad Real 2 3 -0.05790974
## Girona 2 3 -0.04903299
## Guadalajara 2 3 -0.06677319
## Huelva 2 3 -0.06098097
sorted<-sort(km_res_3$cluster)
knitr::kable(sorted, caption="Provincia y cluster")
| x | |
|---|---|
| Barcelona | 1 |
| Madrid | 1 |
| Alicante | 2 |
| Almería | 2 |
| Badajoz | 2 |
| Balears | 2 |
| Castellón | 2 |
| Ceuta | 2 |
| Ciudad Real | 2 |
| Cádiz | 2 |
| Córdoba | 2 |
| Girona | 2 |
| Granada | 2 |
| Guadalajara | 2 |
| Huelva | 2 |
| Jaén | 2 |
| Melilla | 2 |
| Murcia | 2 |
| Málaga | 2 |
| Palmas | 2 |
| Santa Cruz | 2 |
| Sevilla | 2 |
| Tarragona | 2 |
| Toledo | 2 |
| Valencia | 2 |
| Albacete | 3 |
| Álava | 3 |
| Asturias | 3 |
| Bizkaia | 3 |
| Burgos | 3 |
| Cantabria | 3 |
| Coruña | 3 |
| Cuenca | 3 |
| Cáceres | 3 |
| Gipuzkoa | 3 |
| Huesca | 3 |
| León | 3 |
| Lleida | 3 |
| Lugo | 3 |
| Navarra | 3 |
| Ourense | 3 |
| Palencia | 3 |
| Pontevedra | 3 |
| Rioja | 3 |
| Salamanca | 3 |
| Segovia | 3 |
| Soria | 3 |
| Teruel | 3 |
| Valladolid | 3 |
| Zamora | 3 |
| Zaragoza | 3 |
| Ávila | 3 |
También convendría representar algunos boxplot con los valores que poseen las provincias que conforman los clusters en las variables más correlacionadas con las componentes 1 y 2 (al haber tantas en la variable 1 se elegirán Población y PIB) y para la 2, Natalidad y Mortalidad.
data$Grupo<-as.factor(group)
#Boxplot Población
pob<-ggplot(data, aes(x=Grupo, y=Poblacion, fill=Grupo))+
geom_boxplot() +
ylab("")+
ggtitle("Población")
#Boxplot PIB
pib<-ggplot(data, aes(x=Grupo, y=PIB, fill=Grupo))+
geom_boxplot() +
ylab("")+
ggtitle("PIB")
#Boxplot Natalidad
nat<-ggplot(data, aes(x=Grupo, y=Natalidad, fill=Grupo))+
geom_boxplot() +
ylab("")+
ggtitle("Natalidad")
#Boxplot Mortalidad
mor<-ggplot(data, aes(x=Grupo, y=Mortalidad, fill=Grupo))+
geom_boxplot() +
ylab("")+
ggtitle("Mortalidad")
ggarrange(pob, pib, nat, mor + rremove("x.text"),
ncol = 2, nrow = 2)
Con estas visualizaciones se confirman las deducciones realizadas a lo largo de todo el análisis. El grupo 3 que abarca a Madrid y Barcelonaes el cluster con mejor situación socioeconómica, con valores mucho más elevados que el resto en cuanto a PIB, Población y otros indicadores económicos, a parte de tener una buena Natalidad (por haber más población).
Pasando al grupo 2, este incluye provincias no tan excepcionales en términos económmicos como los núcleos del país, pero sí por encima de las del grupo 3. Fundamentalmente son provincias del sur y costa mediterránea junto con muy pocas del interior. Estarían en esta posición por el atractivo turístico ,presencia de buenos puertos para el comercio o proximidad con Madrid Y Barcelona como el caso de Baleares, Málaga, Tarragona, Castellón, Toledo, etc.
Finalmente las provincias del grupo 1 pertenecen a zonas de Castilla La Mancha y Castilla León junto con el norte peninsular, destacando por pobres valores en cuanto a Natalidad, alta mortalidad (por población envejecida y el éxodo rural) así como menor población y PIB debido a la menor presencia industrial y turística.